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1. Introduction 

One of the current problems of interest in aerodynamics is understanding 
and controlling the development of the asymmetric flow which occurs 
when a simple body of revolution is placed at large angle of attack in 
subsonic flow. This phenomenon is of basic scientific interest, since it 
represents a bifurcation from a stable symmetric flow to an unstable 
symmetric (and stable asymmetric) flow as the angle of attack is 
increased. It is also of practical importance in the design of aircraft and 
missiles, since the side forces and yawing moments generated by the 
asymmetrical flow are surprisingly large, and can, in some instances, be 
larger than the normal force and pitching moment acting on the vehicle. 
Designers must take account of these large, nonlinear forces in sizing 
control surfaces and developing control systems. Thus, control of the 
onset of vortex asymmetry can result in improved performance and flight 
safety. 

The phenomenon of vortex asymmetry was first observed experimentally in 
wind-tunnel tests conducted in the 1950‘s, and was then termed 
"phantom-yaw". Early flow-visualization experiments (Ref. 1 ) showed that 
at low angles of attack flow separates from the sides of a circular 
cross-section body and forms a symmetric pair of vortices on the leeward 
side of the body. When the angle of attack is increased beyond a critical 
values, the symmetric vortex pattern changes into a steady asymmetric 
pattern, and large side forces result. In the intervening thirty years a 
great number of experimental and analytical investigations have been 
carried out to document and understand vortex behavior (see, for example, 
Refs. 2-5). It has been found that the onset of vortex asymmetry is 
strongly influenced by several factors, including Reynolds number and the 
state of the boundary layer at the crossflow separation line, model 
smoothness, and tunnel turbulence level. Recent experimental results (Ref. 
6,7), carried out in a low-turbulence wind tunnel, obtained detailed 
surface-pressure distributions over a range of Reynolds numbers ranging 
from fully laminar, through transitional, to “fully" turbulent. These results 
are particularly useful, since they permit delineation of the onset of 
vortex asymmetry from the boundary- layer transition phenomenon. In 
addition, they provide a detailed data base for verification of analytical 
methods. 

Over the past three decades, many theoretical and computational methods 
have been developed and applied in attempting to predict the 
three-dimensional separated asymmetrical vortex flows. These have 
included methods based on the impulsive flow analogy, in which the steady 
three-dimensional flow over a body can be related to the evolution with 
time of the flow over a cylinder in crossflow (Ref. 8-10, for example). 
Also, vortex methods based on conical flow assumptions (Ref. 1 1 ) have 
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been employed. More recently, three-dimensional inviscid vortex-lattice 
methods (Ref. 12), and interactive inviscid/boundary layer method (Ref. 
13) 

have been utilized. However, the strong interaction between the viscous 
crossflow separation and the leeward vortex structure has precluded 
obtaining definitive computational results from these approximate 
methods. As a result, computational methods based on the Navier-Stokes 
equations appear to be required. 

Recently [14], a numerical method based on the thin-layer form of the 
parabolized Navier-Stokes equations [15] was used to compute supersonic 
turbulent flow field surrounding ogive-cylinder and ogive-cylinder- 
boattail bodies at low and moderate angles of attack ( a i 10°). Extensive 
comparisons indicated that the computed results were generally in good 
agreement with experimental measurements [16-181 However, 
discrepancies between the computed and experimental results were seen 
in the regions of experimentally observed crossflow separation. The 
authors of [14] suggested as possible sources of these discrepancies 
between the computed and experimental results the lack of 

circumferential viscous terms within the thin-layer viscous model, and, 
more likely, inadequacies of the algebraic eddy-viscosity model to 
properly treat the regions of separated flow. (Another possible source of 
the discrepancies may have been the marginal computational resolution of 
the leeward-side vortex structures). 

Degani and Schiff [19] improved the PN5 code by adding the cross 
derivative and circumferential viscous terms to the original PNS code and 
modifies the algebraic eddy viscosity turbulence model to take into 
account regions of so called cross-flow separations. The above-mentioned 
sources of discrepancy were investigated. It was found that, given a 
computational grid which provides adequate spatial resolution of the 
leeward-sside vortices, a rational modification of the eddy-viscosity 
turbulence model that is consistent with the physics of the flows extends 
the applicability of the method to flows having large regions of crossflow 
separation. The turbulence model, once modified, was used without further 
changes to compute the flows around an ogive-cylinder body and several 
cones at various angles of attack. The computed results were in uniformly 
good agreement with experimental measurements throughout the flow 
field, for all cases considered, (a s 22.6°). 

As the first part of our investigation we assumed that the flow field is 
conical (but not necessarily symmetric). The validity of the assumption of 
locally conical viscous flow has been demonstrated in Ref. [ 1 4], 
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As a test case for the computations the experiments obtained by Bannik 
and Nebbeling [20] were chosen. In these experiments a cone of 7.5° half 
angle at Mach number 2.94 and Reynolds number 1.372*10 7 was tested up 
to 3 4P angle of attack. At high angle of attack nonconical asymmetric 
leeward side vortex patterns were observed. In the computation, using an 
earlier obtained solution of the above cone for angle of attack of 22.6° and 
at station x=0.5 as a starting solution, the angle of attack was gradually 
increased up to 34° During this procedure the grid was carefully adjusted 
to capture the bow shock and to keep y + near the cone surface smaller 
than 5. A stable, converged symmetric solution was obtained. 

Since the numerical code converged to a symmetric solution which is not 
the physical one, the stability was tested by a random perturbation at each 
point. The possible effect of surface roughness or non perfect body shape 
was investigated too. In all the cases that were investigated the changes 
in the converged solutions were of the order of the surface perturbation. 

At this point of the investigation it was concluded that although the 
assumption of conical viscous flows can be very useful for certain cases, 
it can not be used for the present case and the full marching technique 
should be used (although Peake et al. [21] concluded in their paper that 
conical assymetric vortex pattern is possible). Thus at the second part of 
the investigation an attempt to obtain a marching (in space) solution with 
the PN5 method using the conical solution as an initial data was made, for 
the same cases. 

Since the transition from symmetric to asymmetric pattern probably 
occurs near the tip of the cone the marching should be started very much 
near the tip and a new starting solution should be used instead of the one 
at x = .5 (for conical solutions as described above the location of the 
computed cross section does not matter), and therefore a new conical 
solution near the tip of the cone at x=0.05 is generated for angle of attack 
of 22.6° 

The solution that was obtained using the PN5 method were very much alike 
the solution we got with the assumption of a conical flow field. Thus we 
test the stability of the solution to a random perturbation and ther 
possible effect of the surface roughness as we did with the previous 
solution and again we got a similar behavior. 

All these solutions were done using three different turbulent models. The 
first one is the simple algebraic eddy viscosity model, the second is a 
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modified algebraic eddy viscosity model and the third model was the one 
equation model (or kinetic energy equation model) which uses an additional 
partial differential equation to calculate the eddy viscosity coefficients. 
The model which is used in the present work was formulated by Rubesin 
[22] for compressible flows, based on the Glushko (23) model for 
incompressible flat plate boundary- layer flows. 

For a < 20.6° it was found that the use of the improved eddy viscosity 
model resulted with strong primery symmetric vortices, secondary 
vortices and simple addy viscosity reattachment between them. For the 
case of a = 34° the modified eddy viscosity model was used with the PN5 
code but the solution became unstable. The only solution which shows 
some improvement and was still stable obtained by relaxing the two 
models together. The results were not impressive since the eddy viscosity 
was probably an order of magnitude larger than it should be. The results 
using the one equation model didn't show any improvement as well. 

At this point of the investigation it was concluded that although the PN5 
method code could be very useful for certain cases, it can not be used for 
the present case, of a large angle of attack. The natural step after the use 
of PNS code would be the solution of the full Navier Stokes Equations. 
Since the Navier Stokes Equations are parbolic in time the most common 
method to obtain a steady solution is by marching in time untill 
convergence reached. To overcome the stability problem of the PNS 
solution we had to choose a method that does not suffer from this 
difficulty. The method of flux splitting and upwind spatial differencing for 
the convection terms in the streamwise direction has the advantage of 
having a natural numerical dissipation and better stability properties and 
therefore this method (24) was adopted. 

Using the flux split method we have generated solutions for two test 
cases: The first was the case of the cone that was tested earlier with the 
PNS (but for laminar flow). The second test case for the computation was 
one of the experiments obtained by lemont [7]. In the experiments an 
ogive-cylinder body of 3.5D nose and 7D afterbody. Re = 200,000 based on 
body diameter, and a = 40° was chosen. We assumed that the solution is 
laminar. The solution was made in two steps. As a first step we used a 
grid of 59x50x60 free steam as initial condition and as a second step we 
used a grid of 59x50x 1 20, "where by linear interpolation we used the first 
step as an initial condition. 

The second step was done in three paths. The first one was a time 
marching solution without any disturbance, the second was a solution with 
a perturbation at the nose (small jet perpendicular to the surface) at an 
angle of 90° to the plane of symmetry in the cross section. The third one 
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was a solution where small geometrical disturbance (2% of the diameter) 
was placed instead of the jet and at the same location. 

During the investigation we had to give an answer for two questions. The 
first is, does the solution that we have generated, is the physical one or 
just a methastable solution that will not actually appear. The second 
question is, when we generate a solution using a method of time marching, 
do we reach the seady state solution ? 

Thus a criterion to determine if the flow pattern will actually appear or 
will it be a steady state solution, is needed. We try to develop this 
criterion by using stability analysis. This is based on the assumption that 
solutions which are unstable to small disturbances can not be a steady 
state solution, and metastable solutions appear only for certain initial 
conditions. 

Since the stability analysis is quite complicated and consumes a lot of 
computing time for the three dimension, compressible flow, we needed a 
simple test case where we have several solutions for the same boundary 
conditions. The case that was chosen was Jeffery Hammel flow in a 
diverging channel with large divergence angle [25], (The case of large 
divergence angle, has not been examined for stability yet). Thus the goal of 
this part of the investigation is to develop stability criteria for an 
arbitrary a (where a is the angle between the walls). 

In the stability analysis we limit ourselves to the analysis of linear 
stability of small disturbances such that the nonlinear terms in the 
disturbance equation are neglected. We study the disturbances that fulfill 
the same similarity condition as the base flow, meaning that the 
tangential distributions of velocity of both the distubance and the base 
flow, are independent of the streamwise direction. This assumption is 
equivalent to the case of wave-like disturbance of infinite wave-length in 
the streamwise direction. It was found that the only stable flow pattern is 
the one that an increase of the mass flow rate causes the downstream 
pressure to decrease. 
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2. The General Governing Equations (in body coordinates) 

The general unsteady three-dimensional Navier-Stokes equations, written 
in strong conservation-law form, for Cartesian coordinates can be 
expressed in nondimensional variables as: 


aq aE aF aG 1 aq aR as 
at ax ay az Re ax ay az 


(2.0 


The inviscid flux vectors in Eq. (2. 1 ) are: 


E = E(q) - 


F = F(q) = 


’ pu 


pu 2 + p 



puv 
puw 
(e + p)u 

6 - G(q) - 

pv 


P 

pu 

puv 

pv 2 + p 
pvw 

_ (e + p)v _ 

q - 

pv 

pw 

e 


pw 
puw 
pvw 
pw 2 + p 
(e + p)w 


( 2 . 2 ) 


The internal energy of the gas is defined in terms of the conservative 
variables as 


ej - (e/p) - 0.5(u 2 + v 2 + w 2 ) 


(2.3) 


while the equation of state for a perfect gas with ratio of specific heats 
7 is: 


p/p ■ Cy-Dej » a 2 /? 

The viscous flux trerms in Eq. (2. 1 ) are: 


(2.4) 
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Q - Q(q) - 


■ o ■ 

T xx 

1 

1 

1 

0 

T yx 


0 

T zx 

T xy 

, R - R(q) - 

T yy 

, S - S(q) - 

zy 

i 

M 

J 


V 

Rs . 


T zz 
- Ss . 


(2.5) 


where 


( 2 . 6 ) 


In obtaining Eq. (1.6) the Stokes hypothesis was used: x = -2ji/3. In Eqs. 
(2.1M2.6) the Cartesian velocity components u ,v ,w are made 
nondimensional with respect to a (the freestream speed of sound), the 

density p is normalized by p w and total energy e is referenced to 

a 2 
r °o oo 

We introduce a general transforamtion as follows: 


L 


I - £(x,y,z,t) 
Ti(x,y,z,t) 
c = £(x,y,z,t) 


(2.7) 
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Subject to the transformation, Eq. (2.1) can still be expressed in a strong 
conservation-law form as: 


A A A A AAA 

aq 3E 3F 36 1 6Q 3R 3S ^ 

it* IT 3t 1 K Re 3 t 3 n K 

(2.8) 

where 


t 

i 

A 

q - q/J 


E - (5„E' * F ♦ ? Z G- ♦ f t q)/J 


A 

F = (t] x E + T)y F + r| z G + T^q)/J 


6 = + Cy R + + 

i 

i 

1 

(2.9) 

and 


Q - Qf (q,q|) + QH (q,q n ) + fc(q,q,0 - (f x Q + f yR + £ Z S)/J 

(2.10) 

R = R^(q,q|) + RH (q,^ ) + &(q,q{) = (n x Q + *ly R + flz s > /J 

(2.11) 

S « S^(q,q|) + S 1 ! (q,q n ) + sC(q,q^) = ({ X Q + < y R + £ Z S)/J 

(2.12) 


The Jacobian of the transformation, which appears in Eqs. (2.8M2. 12), is 
defined as: 


J ' 1 = ^ * «l (z n *( ' Y{ 5 * z 5 (x n ' Vi’ 


(2.) 3) 
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3. The P.N.5. Code 


The P.N.5. code deals with the steady state solution. Thus the variables 
are independent of time. The parabolized Navier-Stokes equations are 
obtained from Eq. (2.1) by a generalized independent spatial variables that 
map the physical x,y,z space surrounding a body into a rectangular ^ ,q, {; 
computational region: 

| - % (x) - streamwise (marching) coordinate, 

r\ = q(x,y,z) - spanwise or circumferential coordinate, 

£ = £(x,y,z) = normal coordinate. 


which maps the body surface into the £ - 0 plane, and by neglecting all 
streamwise derivative, 3/3£, within the viscous terms, and by modifying 
the streamwise flux vector to permit stable time-like marching of the 
equations downstream from initial data. Following [15], we introduce the 
subsonic sublayer approximation, and the resulting parabolized 
Navier-Stokes equations can be written as: 


AAA 

3E S 3F 3G 
3f + 3q + 3{ 



3 a a 3 a a 

— (Rn + RC) + — (S^ + sC)] 

3q 3^ 


( 3 . 1 ) 


The modified streamwise flux vector in Eq. (2. 1 ) is 


A 

E s - (| X /J)E S - (f x /J) 


pu 

PU 2 + Pg 

puv 

puw 

(e+p s )u 


( 3 . 2 ) 


where p s = p for supersonic flow, and p s is defined from 3p / 3£ = 0 for 

subsonic flow in the viscous layer adjacent to the body surface. By 
evaluating p s in this manner, Eq. (3.1) can be stably marched in the X 

direction for all flows where u > 0; that is, for flows without streamwise 
reversal (see [15] for associated stability analysis). 



The viscous flux vectors in Eq. (2.1) R 7 ! , R$ , S 7 ! and S 7 ! are given in 

the appendix. 


3.1 Numerical Algorithm 

The numerical algorithm used to march Eq. (3.1) downstream is an 
implicit, non-iterative, approximately factored finite-difference scheme, 
which is analogous to the one developed by Beam and Warming for the 
solution of the unsteady Navier-Stokes equations. The marching algorithm 
is derived in the same manner used by Schiff and Steger [15], but the 
viscous cross-derivative term 3R ^/dq dS 7 ) /a$ in Eq. (3.1) cannot be 
treated implicitly, and, following Beam and Warming [26], are evaluated 
explicitly. The resulting algorithm can be written in so-called delta form 
as 


(A S J + (l-cOAtf^ BJ - Re’ 1 ^ N^KAgl)" 1 - 
* [A S J + (l-oOA^cJ - Re^ffojAqJ 

-(Agi-Agi-^qj + attgJ-EgJ’ 1 ) - (1-a)A|{$ n fo x j+1 (E/J)i + 


+ q y i+ 1 (F/J)J + ii 2 j + 1 (S/J)j] + a^ x j + 1 (E/J)j + C y i +, (F/J)j + 

+ C 2 j+1 (G/J) j ] - Re" 1 [6^ + 6^ (rC)J + 6^)1 + 

♦ tySty] - BRe^lfyARty " 1 + 1 ]( - [({ X /J)j +1 E p J 

- <* X /J)J E p i - 1 ] ♦ Dqj ♦ 0 (Ae) 1+3a ( 3 . 3 ) 


In Eq. (3.3) a is set equat to 0 for first-order accuracy (Euler implicit 
method) and a = 1/3 for second-order accuracy (3-point backward 
differencing). Similarly, in the viscous terms e = 0 for first-order 
accuracy and 8 = 1 £or second-order accuracy. The Jacobian matrices of 
the flux vectors, A s , B, and C are obtained from local linearization [28] of 

A /s/ 

E s , F and G. The Jacobian matrix li is obtained from local linearization 
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A L ^ 

[28] of 5^ and, in an analogous manner, the Jacobian matrix N is obtained 
from linearization of R 7 ! . The symbol ~ indicates that the matrices are 
evaluated using flow variables q located at jA£ and metric quantities 
at (j + 1 )A£. The term Dqj is a fourth-order explicit dissipation term, 
defined as: 


DqJ - e e AJ s k i ( 1 /J) J l(V n A^) 2 (Jq)l + (V^A^) 2 (Jq) j] (3.4) 


which is added to the algorithm to suppress high-frequency oscillations. 
Linear stability analysis indicates that e e in Eq. (3.4) must be less than 

1/16 to ensure stability of the algorithm. Although adding implicit 
smoothing terms within the operators on the left-hand side of Eq. (3.3) 
overcomes the linear stability limit and permits the use of larger values 
of £ e , no such implicit smoothing was used. While the use of implicit 

smoothing tends to stabilize the numerical method, the added smoothing 
terms can be larger than the viscous terms of interest, and thus can 
degrade the accuracy of the solution. 

Equation (3.3) contains all viscous terms and viscous cross terms 
applicable to the parabolized Navier-Stokes equations. The thin-layer 
viscous model form of the equations, previously used by Schiff and Steger 
[15] can be obtained by neglecting all viscous terms except those solely in 
the normal $ direction. Thus, by dropping the terms ft, R, S 7 ] , AR S , and 
AS 7 ! from Eq. (3.3), the thin-layer algorithm can be obtained. 
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4. The Conical Solution using the PN5 Method 

In general, the initial data for the marching method must be supplied 
from an auxiliary computation. However, when treating the flow over 
conical or pointed bodies this is not necessary. As outlined in [15], for 
inviscid flows about conical bodies a conical grid is selected and the flow 
variables are initially set to free-stream values. The solution is marched 
one step downstream from an initial station, and the resulting flow 
variables are then scaled to place the solution back at the original station. 
The process is repeated until no change in the variables is observed with 
further marching. The flow variables are then constant along rays of the 
flow field, and a conical solution has been generated. Upon assuming flow 
variables within the viscous layer to also be constant along rays, the same 
procedure can be used to generate viscous conical solutions. Although 
viscous flow cannot be strictly conical, the assumption that the flow 
variables are locally conical along rays is reasonable when treating 

high-Reynolds-number flows. The validity of the assumption of locally 
conical viscous flows has been demonstrated in [14]. The 
marching-stepback procedure was utilized to generate the conical 
solutions in the present work. 


i 
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5. The Flux Split Scheme Applied to the Thin Laver Approximation. 


The flux split method is a time marching method and thus the Eq. to be 
solved is Eq. (2.8), The thin-layer approximation means to neglect the 
viscous terms in the streamwise and the spanwise direction and keep only 
the viscous terms in the normal direction. Thus the Equation to be solved 

is: 


A A * £ * 

aq 3E 3F 3G 3s. 

— + — + — + 77 - Re' 1 7T 
3t 3| 3q K ^ 


(5.1) 


where 5 now is: 


S 



’ ^<x 2 * iy 2 * <z 2) \ * t 1 ' 3 «x u { * 

VLCCx 2 + Cy 2 * ^ 2 2)v { * l 1 ^ 3 <tx u { * <<j\ * < 2 w {Xy 

* H <{ x 2 * <y 2 * tz 2)w { * ^ /3 «x u { * V? * 

((^ x 2 ♦ (J * ^ 2 j )[o.5h(u 2 + v 2 + *> { ♦ xPr" 1 (y- 1 )' ’ (a 2 )^] ♦ 
* li/3 ({ x u ♦ { y v - { 2 w)({ x u { ♦ ♦ C 2 w f )) 


(5.2) 


5. 1 Numerical Algorithm 


The Split Flux method was developed to overcome the difficulty of 
stability. This goal could be achieved by combining two methods: The first 
is by using two steps factored scheme which was proved by Ying [28] to be 
unconditionally stable (differently from the three factored schemes), and 
the second is by splitting the Flux in £ direction according to its 
eigenvalues. It was proved by Ying [28] that the stability A of the scheme 
depend on the eigenvalues of the matrix E. The splitting of E is equivalent 
to adding a viscosity that stabilize the scheme. 


The two steps implicit factored scheme with splitting of the E flux 
for the thin layer approximation using the upwind approximation in the ^ 
direction and central differencing in the t \ and $ direction can be written 
as follows: 



1 5 


[I * h6 ? b (A*) n ♦ h6^C n - hRe' ' 1 M n J - 0,^1 

. II ♦ h5 t f (A") n * h« n B n - D ttn 1 iq n - 
- -At(5| 6 t(E*) n - E m + 1 ♦ 5| f KE") n - E„'l * 6 r| ^ * 

* 8^(3" - GJ - Re- 1 6 { (S"-SJ} - Dglq" - qj 


where h » eAt. 8, 8j= b » 8 ^ are central, backward and forward three points 

difference operator A 1 , B, C, M are obtained from the local linearization 
of E F, 6, and S. E* is the flux E which has been split into E + and E' 
according to its eigenvalue in the following manner: 


A 


E = 


8i(Xi,A4,A5) 


( 5 . 4 ) 


and 

2(y-1)xr + X4” + xs 1 
2(y-i)xi ± u * X4*u + xrir 
E 1 - p/2y = 2(y- 1 )xrv + xA + x 5 4 v' 

2(y-l)xrw + x^w + X 5 4 W" 

24. * X4*/2 l(U*) 2 * (V‘) 2 ♦ (w*) 2 ] * M ± [(U') 2 ♦ (V) 2 * 

♦ (W") 2 1 * C m (5.5) 


C III 


(3-y)(X4 ♦ xs)c 2 

2(7-1) 



and D v and D e are the numerical dissipation terms which are given as 

combination of second and fourth differences. The smothing terms have the 
following form: 


D e(i] " (AtXJ~ 1 {£2^p(B)p6 + 846 3 3 } 1 j 

1 + B 


1 + p 


d i|t| 58 (At)j " 1 {e 2 6 p(B)p 6 + 2.5 £46 



where 


(5.6) 


P 


1S 2 P| 

IO* 6 2 )p| 


(5.7) 


A 


A 


P is the nondimensional pressure and p(B) is the spectral radius of B. 


6. . Turbulence Models. 

6. 1 The Eddv-Viscosity Model 

In this discussion we adopt the dimensional notation of Baldwin and 
Lomax [29], The resulting coefficients can be nondimensionalized for use 
in Eq. (3.3) by normalizing them by their free-stream (laminar) values. 

For laminar flow computaitons the coefficient of molecular viscosity 
ji ■ ji 2 * s obtained from Sutherland’s law and the coefficient of thermal 

conductivity k is specified, assuming a constant Prandtl number, as 
K 2 /Cp = jiL/Pr. For turbulent-flow computations the laminar flow 
coefficients are replaced by 


M- ■ h* H 


K/Cp = jJL^/Pr + (jL t /Pr t 


( 6 . 1 ) 



17 


The turbulent viscosity coefficient is computed using the isotropic, 

two-layer, Cebeci-type, algebraic eddy-viscosity model reported by 
Baldwin and Lomax [29], 


In the Baldwin-Lomax formulation ji^ is given by 

lit = { ^ inner ' y i yc 
Wouter » y > tic 


( 6 . 2 ) 


where y is the local distance measured normal to the body surface and y c 

is the smallest value of y at which the values from the inner and outer 
region formulas are equal. Within the inner region 


<IH >tnner " Pt*M 

where 


i = ky[l 


e-ty/Mj 


Itol is the magnitude of the local vorticity vector, and 


(6.3) 


(6.4) 


y - 


(6.3) 


In the outer region, for attached boundary layers the turbulent viscosity 
coefficient is given by 


Pouter ^vp P^weke ^kleb^ 


(6.6) 


In Eq. (6.6) K and C cp are constants, F k | eb is the Klebanoff intermittency 
factor, and 
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F wake yrnax F niax' 

where F max is the maximum value that the function F(y), defined as 


(6.7) 


-(u + /A + ) 

F(y) = |co| y[1 - e ] (6 - 8) 

takes in a local profile, and y max is the value of y at which F max occurs. 

The constants appearing in Eqs. (6.1H6.8) were determined in [29] by 
requiring the boundary- layer profiles computed wilth the model to be in 
agreement with those determined using the Cebeci [30] formulation. The 
values were determined to be 


Pr = 0.72, 

O 

tl 


Pr t = 0.9, 

K = 0.01681, 

(6.9) 

A + = 26, 

Cep-' -6- 



6.2 Modified Addv Viscosity Model 

The major difficulty encountered in applying the Baldwin-Lomax 
turbulence model to bodies with crossflow separation is that of properly 
evaluating the scale length y max and in turn, of determining (ji t ) outer for 

boundary- layer profiles in the crossflow separation region. This difficulty 
becomes apparent upon considering the behavior of the function 
F(y)[Eq.(6.8)] along two rays, one located on the windward side at $ = $ j 




Fig. 6.1 Behavior of F(y) at large incidence, (a)* = +, (windward side), 
(b) ♦ = * 2 (leeward side). 
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and the other on the leeward side at + = * 2 . The functions are shown 

schematically in Figs. 6.1(a) and (b), respectively. On the windward side 
the attached boundary layer gives rise to a profile of F(y) which has a 
single, well-defined, peak, as shown in Fig. 6(a). Thus, the determination 
of F max^ *1 ymax^P’ and F wake * i ) Is straightforward. However, on the 

leeward-side ray [Fig. 6.1(b)], in addition to a local peak in F(y) in the 
attached boundary layer at y 2 - a, the overlying vortex structure causes a 

larger peak in F(y) at y 2 - b. As originally implemented, the computer code 

searches outward along each ray to determine the maximum in F(y), and 
would, in this instance, select the peak in F(y 2 ) occurring at y 2 = b. The 

choice of the peak at y 2 - b results in a value of F wake ( * 2 ) and, in turn, a 

value of the outer layer eddy-viscosity coefficient (Router wh ’ ch is 

much too high. The resulting value is at least one order of magnitude 
larger and can be as much as two orders of magnitude larger than the value 
of (MtWer resulting from evaluating F wake ( * 2 ) from the peak at y 2 = a. 

Thus, in general, the computed eddy-viscosity coefficient in the crossflow 
separation region behind the primary separation point will be too high. 
This will cause the details of the computed flow to be distorted or washed 
out. In particular, the primary vortices will be smaller than those 
observed experimentally and the primary separation point will be located 
closer to the leeward symmetry plane. In addition, the secondary 
separation and secondary vortices will not appear in the computed flow. 

To eliminate these difficulties we have modified out implementation 
of the turbulence model. At each axial station the code searches radially 
along successive rays, sweeping from the windward to the leeward plane 
of symmetry. Along each ray the code sweeps outward to find the first 
peak in F(y), and cuts off the search when the peak is reached. To prevent 
the selection of extraneous peaks which might be caused by a nonsmooth 
behavior in F(y), a peak is considered to have been found when the value of 
F(y) drops to 90 % of the local maximum value. Choice of F max in this 

manner will exclude the second, spurious, maximum [see Fig. 6.1(b)], 

For most rays in the crossflow separation region the two peaks in F(y) 
are spaced far enough apart that the logic described above will select the 
first peak. However, this is not true for rays in the vicinity of the primary 
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separation point (and to a lesser extent for rays in the immediate vicinity 
of the secondary separation point). Along these rays the vortex feeding 
sheets lie close to the outer edge of the attached boundary layers: the 
peaks in F(y) merge. Under these conditions the code would choose a value 
of y max near the top edge of the feeding sheet. Consequently, a further 

test is applied. On each ray (except the ray on the windward plane of 
symmetry) a cutoff distance is specified in terms of y max from the 

previous ray, i.e., y cu t 0 ffW = c Ymax (♦"*♦)* where c is a constant 
chosen equal to 1.5. If no peak in F(y) is found along a ray for y^y cu toff the 
values of F max and y max are taken as those found on the previous ray. In 

this manner a physically reasonable value of the eddy-viscosity 
coefficient will be chosen for those rays close to the crossflow 
separation points. 

It Is readily apparent that conditions withing the boundary layers 
which leave the body at the primary separation points are related to the 
conditions within the boundary layers on the windward side of the body. 
Further, it is physically reasonable to expect that the boundary- layer 
quantities vary smoothly circumferentially around the body. Thus, 
specitying the cutoff distance in terms of the values on the previous ray, 
and taking the values of y max and F max from those of the adjacent ray, 

allows the model to be applicable in a rational manner over a wide range 
of local flow conditions, and in particular, for varying local Reynolds 
numbers. 

The various coefficients appearing in the Baldw in-Lomax turbulence 
model were varied to assess their effect on the computed boundary- layer 
profiles. The best match with experimental measurements was obtained 
with the coefficients set at the values suggested in Eq. (6.9). 
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6.3 One Equation Model. 

This model has one turbulence source function H^: 
H k = ^S2 - 2pkD/3 - Cn k k/1_2 

where 

S 2 = (u y +v x ) 2 + 2[u x 2 +v y 2 + (v/r) 2 ] - 2D 2 /3 
D = u x + v y + v/r 
- ccji^ R t H (R t /Ro) 

R^ » p V'R L/m 


H(R) 


4 


R 

R - (R-0.75) 2 
1 


for R < 0.75 
for 0.75 <R< 1.25 
for 1.25 < R 


a - 0.2 , C - 3.93, Rb - 1 10, P rk - 2.5 


Here k is the turbulence kinetic energy, k = pU^' 
length scale as specified by Glushko [23]: 


y/5 for 0 < y/5 < 0.23 

L/5 = (y/5+0. 37/2.61 for 0.23 < y/5 < 0.57 

(1.48-y/5)/2.52 for 0.57 < y/5 < 1.48 


(6.10) 


(6.11) 


‘j 1r /2p, and L is the 


(6.12) 
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The total kinetic energy diffusivity ji^ is given by: 

t*k " m * ViPr 

k 


The total viscosity |X and the thermal energy diffusivity are 


H - V-x + H 

K/Cp - /Pr + |it/Pr t 


( 6 . 13 ) 


( 6 . 14 ) 
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7. The Test Case for the Stability Criterion. 

The test case that was chosen is the Jeffery Hammel flow in a divergent 
channel since for that case we have an analytic solution which provides 
two different patterns of flow for the same boundary conditions. 


7.1 The flow field 



Fig. 7.1 - Configuration of the problem. 


Statement of the problem. 

We study a divergent channel, infinite in the z direction (Fig. 7. 1 ) in 
cylindrical coordildnates, r, 8, z. r=0 is the (virtual) intersection of the 
two walls. A Newtonian, viscous incompressible fluid flows in the channel 
as a result of inertia forces or an external pressure gradient which is a 
function of r and 8. We assume that the flow is isothermal and has only a 
radial component V r Thus 


v z = V 0 = 0 everywhere in the flow 


(7.1) 


where V z and V 0 are the velocity components in the 8 and z directions 
respectively. 

The boundary conditions are the no-slip conditions: 



at the walls 


(7.2) 


V 


r 


= 0 


where Vr is the velocity component in the r direction. 

We consider the mass flow rate or the pressure gradient at the walls 
as given (compatibility conditon). 

By defining nondimensional velocity and time as: 


u = V /u 
u v r u max 


T = t 


Umax 


(7.3) 

(7.4) 


we obtain a Re number: 


Re = u max ry,v 


(7.5) 


From the incompressible Navier-Stokes equation we get the time 
dependent equation for the flow field: 


Re (3U/at) - (3 2 U/30 2 ) + 4U + ReU 2 + K 


(7.6) 


where 


3 2 U 

K - K(t) - -c(t)/Re - - — t 


walls 


(7.7) 
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The pressure equation is 

P - P o (0,t) + 2p ReU + KRe <pv2/2r«) 

and the pressure gradient along the walls (ap/dr) w , is: 

(ap/ar) w - (-K Re) pv/r 

The mass flow rate M, for a unit depth of the channel is: 

M-Re^Jude (7.10) 

cc 

where ji is the dynamic viscosity. 

7.2 Solution for the Steady State Flow Field. 

The solution for this case was presented by Shapira, Degani and Weiss 
[25] as follows: 


(7.8) 


(7.9) 


U’ = ± / (2Re/3) ( 1 -U) (p 2 - U) (p3 ‘ U) (7.11) 

where 


p2 = a + d 


p3 = a-d 


(7.12) 
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and 


a - -( 3 /Re + 1 / 2 ) 


d = 1 /( 2 Re) V 3(12 - Re( 4 k + Re + 4 )] 


(7.13) 


7.2. 1 The steady state solution for -(Re+6)/3 < k < - Re/4 + 3/Re - 1 is: 
The solution in the above range is given by: 


U ■ 1 - (1~p2) Sn 2 (9 /m , ki) 


( 7 . 14 ) 


where 


1 -p2 6 

ki 2 = M 2 = 


1 - p3 


Re( 1-p3> 


and Sn(0/X 1f k 1 ) is the Jacobian elliptic function of the first kind, and 


I 


b-m| 


% 


o (1-ki 2 sin 2 |) 1/2 


M F(e,k,) 


( 7 . 15 ) 


For the channel angle a and the corresponding mass flow rate we can 
get a number of solutions which fulfill one of the following rules: 
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a = 2(n ai + (n-1) « 2 ) ; M = 2(nMi + (n-1) M 2 ) 

or 

oc * 2(n <x\ + (n 0 C 2 ; M»2(nMi + nM 2 ) 

or 

oc = 2(n ai + (n+1) 0 C 2 ) ; M = 2(n Mi + (n+1) M 2 ) (7.16) 

where n- 1,2,..., and a < ti . 

and 


0(1 o (1-ki 2 sin 2 f) l/2 


M F(tk,) 


(7.17) 


where F(£i, ki) is the incomplete elliptic function of the first kind. The 
mass flow rate passing through this element is: 


Rep. 

Mi = —.[(k^ - 1 + p2)cci + X 1 O-P 2 ) E(£i, ki)J (7.18) 

K I 

where E(f \, ki) is the incomplete elliptic function of the second kind. 
The angle of the element 0 > U > P 2 is: 


n/2 of 

( 1 -ki 2 sin 2 $) l/2 “ ! ' ,F(k|) '“ l 


(7.19) 


where F(k pis the complete elliptic function of the first kind. The mass 
flow rate passing through this element is: 


Reii 

«2-M,- I(k2,- 1 +p 2 )F(k 1 ) + (1-p 2 )E(k 1 )] 

ki 2 


(7.20) 
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where E(kj) is the complete elliptic function of the second kind, and 

■ arcsln (1/ l/u-p 2 ) ) (7.21) 


7,2.2 The steady state solution for k > - Re/4 ♦ 3/Re + 1 
The solution in the above range is given by : 


U - 1 - L( 1-Cn (6/*2, k 2 )) / ( 1 + Cn (8/ A2 , k 2 )) • (7.22) 

where 

L2 = 3(K + Re + 4) / Re (7.23) 

and 


X2 2 = 3 / (2 Re L) 

k 2 2 = (L + 3/Re ♦ 3/2) / (2 L) (7.24) 

and Cn(8/A2,k 2 ) in the Jacobian elliptic function of the second kind 

0 * x 2 f 

o (1-k^ sin 2 1) ,/2 


N2 F<e,k 2 ) 


(7.25) 
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The half of the channel angle is: 


a/2 - \2 I 


%2 




o ( 1 -k 2 2 sin 2 £) u2 


*2 F(^2,l<2) 


(7.26) 


The mass flow rate is: 

M = Re|i [( 1 -L)a + 4Lx 2 (E(fc,ki) - L2)] (7.27) 

where 

L2 = / L(L-4L k^ + 2L + 1 ) / [L(L+ 1 )] (7.28) 

and 

%2 = arcos ((L- 1 ) / (L+ 1 )) (7.29) 


7,3 The stability equation. 

Next we examine the linear stablility of the flow patterns subject to 
perturbations that vary as 1/r and are functions of 6 and t. 

Let us define a nondimensional solution: 

U(8,t) = uo(0) + w(0,t) (7.30) 


which was obtained by multiplying the solutions of the previous section by 
r, where Uo is the base flow that is analyzed for stability, and W is the 
perturbation. 

Substituting Eq. (7.30) in Eq. (7.6), subtracting eq. (7.6) (written for 
Uo) and neglecting the nonlinear terms yields the following linear equation 
for small disturbances: 


aw/at = (vr + 4W) / Re + 2UqW 


(7.31) 
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Let us define a disturbance of the type: 

W(B,t) = W( 8 ) exp( p t) (7.32) 

From Eq. (6.3) we obtain: 

(4/Re ♦ 2Uo)W + WVRe - W = 0 (7.33) 

With the following boundary conditions: 

W = 0 at the walls (7.34) 

Eq. (7.33) and the boundary conditions (7.34) are a system of 
eigenvalues p and eigenfunctions W(8). 

The base flow pattern is stable if: 



Let us define the following approximation: 

tyi-1 + Wn+1 " 2W n 

Wn- n * 0, 

h 2 


(7.36) 


where h is the difference between two discrete points and Wn is the value 
of W at 8 = n*h. 
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Applying Eq. (7.36) to Eq. (7.33) yields: 

sW * Qn (Uo) Wn + *ib Wnt ' ’ p Wn ■ 0 

where 


Q n (Uo) - (4/Re + 2Uo 


- 2/h 2 Re ) 


8 = nh 


and 




= 0 


Utilizing Eq. (7.37) enable us to turn the eigenvalue problem into: 

det( A - pi) = 0 
where I is the unit matrix, 
and 


A 


Qt 1/h 2 Re 

0 

1/h 2 Re (fe 


0 


Q 1/h 2 Re 
Vi 


(7.37) 


(7.38) 


(7.39) 


(7.40) 


1/h 2 Re Q n 


(7.41) 
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8. Results and Discussion 


8.1 The Stability Analysis 

8.1.1 Decreasing pressure in the downstream direction. 

In the case where the pressure decrease is in the flow direction two 
counteracting mechanisms exist: The first is the pressure gradient that 
tries to accelerate the flow. The second is the deceleration and 
accompanying pressure rise due to the radial spreading of the flow. 

In Fig. 1 contours of constant ct are plotted in the plane of M (the 
mass flow rate) vs. (ap/dr) w (the pressure gradient on the walls). As a 

result of the two opposing mechanisms mentioned above there is a critical 
value of the mass flow rate M c for each a that behaves in the following 

manner: For M < M c an increase in the mass flow rate, will cause the 
downstream pressure to decrease. When M > M c an increase in the mass 
flow rate, will cause the downstream to increase. 

Another result that can be deduced from Fig. 1 is: For any pressure 
gradient on the walls (ap/dr) w , we can find two values of mass flow 

rate M. One is greater than M c and one is smaller. 

Curves of M c are depicted In Fig. 1 and Fig. 4. The accompanying K and 
(dp/dr) w are given in Fig. 6 and Fig. 3. 

The stability method, described in Paragraph 7, leads us to the 
following conclusion: When M > M c the flow is unstable, and when M < M c 

the flow is stable, but when M > M c , the mechanism that limitss and 

stabilizes the flow, is missing. 

We also see from Fig. 1 that for any a the pressure gradient and the 
mass flow rate M are bounded. In Fig. 3, we see the line that separates the 
zone of two solutions and the zone where no solution of the type Vg = 0 

exists. Fig. 2 depicts graphs of pressure gradients vs. a for constant Re. 
We see that each curve ascends to the contour depicted in Fig. 3 becomes 
tangent to it, and then returns back into the two solutions range. From the 
stability analysis we found that the solution is unstable in the range of 
ascent and stable afterwards. 
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In Fig. 4, zones of existence of the solutions and their stability are 
depicted. For any a there is a stable zone of low mass flow rate, an 
unstable zone of higher mass flow rate, and a zone where no solution of 
the type V = 0 exists. 

In Fig. 5 contours of constant a are ploted fora between 5° and 60°. 
From Fig. 5 we see that the contours of constant a are slightly convex. 
Thus for any given a it is possible to find two constants which 
approximate the relation between K and M. In order to approximate the 
velocity profile, using these three constants, the Re number is needed. In 
Fig. 6 contours of constant Re in the a vs. K plane are given. From Fig. 6 and 
Eq. (7.14) or Eq. (7.22) the pressure gradient and the velocity profile can be 
calculated. 

In Fig. 6 the stability line in the K vs. a plane is depicted, so that 
below this line the solution is unstable, and stable above it. 


8. 1 .2 Pressure gradient in the upstream direction. 

In the case where the pressure decreases in the upstream direction 
three patterns of flow can be defined (as expresses in Eq. (7.16)): The first 
is a syummetric pattern of flow without inverse flow near the walls 
(positive shear stress on both walls). The second is an asymmetric pattern 
of flow with an inverse flow near one wall (positive shear stress on one 
wall and negative shear stress on the other). The third is a symmetric 
flow pattern with an inverse flow near both walls (negative shear stress 
on the two walls). 

It is obvious that for any pattern we can find several zones of 
negative and positive velocity that obey the rules appearing in Eq. (7.16). 
As mentioned in Paragraph 7 when K = Kp, the range of inverse flow is 

reduced and the three patterns of flow become the same. 

In Fig. 7 curves of constant a are plotted in the plane of M (mass 
flow rate) vs. (dp/dr) w (the pressure gradient on the walls), for the case 

where only one range of positive velocity exists. In this case, for anya we 
have three branches with a common point at K = K p . The left branch is for 

the case without inverse flow, the centre branch for inverse flow near one 
wall and the right branch describes the case with inverse flow near both 
walls, thme left branches are the continuation of the curves obtained when 
the pressure decreases in the downstream direction. 
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From Fig. 7 we conclude that the contours are qualitatively similar 
for the full range of a. Fig. 8 shows curves of a = 20° in the plane of 
(dp/dr) w vs. M, for single double and triple positive veleocity zones. 

From Fig. 7 and Fig. 8 we also conclude that in the case of the 
velocity profile without an inverse flow (left branches in Fig. 7) the 
solution is unique. In patterns of flows where zones of negative velocity 
appear, for every pressure gradient we have at least two solutions for the 
mass flow rate. Fig. 9 describes the zones where one solution and multiple 
solutions exist, in the plane of a vs. (dp/dr) w . 

In Fig. 7, we see that for pressure decreasing in the upstream 
direction we can have higher mass flow rate than we had in the case of 
pressure decreasing in the downstream direction. From Fig. 7 and Fig. 8 we 
understand that it is possible to define the ranges of mass flow rate and 
the range of pressure gradient for each flow pattern. Beyond these range 
this pattern cannot exist (similar to Fig. 9). 

From the stability analysis we found that all the patterns of flow 
where the pressure gradient is in the upstream direction are unstable. We 
found that there are cases where the pressure gradient is in the upstream 
direction and increasing the pressure gradient causes an increse in the 
mass flow rate. Analysis based on physical reasoning could lead us to the 
conclusion that these patterns are stable. The stability method, described 
in Paragraph 7, shows however that the flow is unstable. This instability 
can be explained by the combinations of unstable elements. The instability 
expresses itself in the fact that, in such elements a decrease in the mass 
flow rate causes an increase in the pressure gradient. 



8 2 Conical Solution 


The Schiff-Steger PN5 code [aalhas been modified to allow computation of 
conical flowfields around cones at high incidence. The improved algorithm 
of Degani and Schiff [aa] has been incorporated with the PN5 code. This 
algorithm adds the cross derivative and circumferential viscous terms to 
the original PN5 code and modifies the algebraic eddy viscosity turbulence 
model to take into account regions of so call cross-flow separation. 
Assuming the flowfield is conical (but not necssarily symmetric) a march- 
ing stepback procedure is used: the solution is marched one step down- 
stream using improved PN5 code and the flow variables are then scaled to 
place the solution back to the original station. The process is repeated 
until no change in the flow variables is observed with further marching. 
The flow variables are then constant along rays of the flowfield. The 
experiments obtained by Bannik and Nebbeling [20] were chosen as a test 
case. In this experiments a cone of 7.5° half angle at Mach number 2.94 and 
Reynolds number 1.372*10 7 was tested up 34° angle of attack. At High 
angle of attack nonconical asymmetric leeward side vortex patterns were 
observed. In the first set of computations, using an earlier obtained 
solution of the above cone for angle of attack of 22.6° and at station x=0.5 
as a starting solution, the angle of attack was gradually increases up to 
34°. During this procedure the grid was earful ly adjusted to capture the 
bow shock and to keep y+ near the cone surface smaller then 5. A Stable, 
converged symmetric solution was obtained. It was the first time 
that a numerical solution was obtained for an angle-of-attack-to- 
half-cone-angle ratio as high as this case. 

Since the numerical code is perfectly symmetric, in the second set of 
computations a random perturbation of about 1% of the local flowfield 
variables was introduced at each point of the flowfield and the marching 
stepback procedure was continued till a new converged solution was 
obtained. It was found that in all cases tested the solution converged back 
to the original symmetric solution. By increasing the initial perturbation 
above 2% - 5% of the local flowfield variables the numerical solution 
became unstable and did not converge. 

In the third set of numerical experiments the possible effect of surface 
roughness or non perfect body shape was investigated. The cross section 
of the cone has been changed randomly up to 1% of the original diameter 
and kept so until a new converged solution was obtained using the 
marching stepback procedure. It was found that in all cases investigated 
the changes in the converged soltuions were of the order of the surface 
perturbation. 
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8 3 Marching Solution 

8.3. 1 Algebraic T urbulence Model 

As mentioned earlier, Bannik and Nebbeling [20] found in their experiments 
that at high angle of attack a nonconical asymmetric vortex pattern is 
obtained. Therefore, for these cases the full marching technique should be 
used. Moreover, since the transition from symmetric to asymmetric 
pattern probably occurs near the tip of the cone, a new statring solution 
should be generated instead of the one at x=. 5 (which was used for the 
conical solutions). A new conical solution near tip of the cone at x=0.05 
was generated for angle of attack of 22.6° Figs. 10-11 show results 
with the old algebraic turbulence model and Figs. 12-13 show results 
using the modified algebraic model. From the velocity vector plots one 
can see that the use of the modified model resulted with strong primery 
symetric vortices and secondary ones. The results of the old model show 
just one weak pair of primery vortices. Using these results another 
conical solution for angle of attack of 3^ at x=0.05 was generated. Figs. 
14-15 show the results for this case using the old turbulence model At 
this point the modified algebraic model was introduced but the solution 
became unstable. The only solution which showed some improvement in 
comparison to the solution obtained with the old turbulence model and 
was still stable obtained by relaxing the two models together. A solution 
which obtained using a ratio of 90% of the value of eddy viscosity from 
the improved model and 10% of from the old model is shown in Figs. 16-17. 
The results are not impressive and one has to remember that although 
the ratio is 9:1 the eddy viscosity is probably still too high. 


8.3.2 One Equation Turbulence Model 

For flowfield where the mean flow changes so rapidly that the turbulence 
cannot remain in equilibrium with mean motion, an algebraic model might 
not be sufficient. In cases where the region of separation is relatively 
small compare to the upstream boundary layer thickness a one equation 
turbulence model can be used (see paragraph 6.3). Although the use of this 
model improved the simulation of separated flows in some cases (see for 
example Ref. [31]), it did not improve the unstable behaviour of our 34° 
angle of attack case. 


8.8.3 Asymmetric solution for cc =34° 


Using the combined old and modified turbulence models , as mentioned 
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above, an asymmetric solution was obtained by perturbing the initial 
flowfield. Starting with a conical (symmetric) solution near the tip of 
the body, the PNS code was used to march down the body at 34° angle of 
attack. During the first ten steps a small perturbation of about 1% of the 
local value of the symmetric solution was added in random. Figure 18 
shows the density contours at different stations down the body and Figure 
19 shows the velocity vectors for the same stations. 

The results shows large asymmetry of the solution (And it is much larger 
than the initial perturbation) and the primery vortices are not straight or 
parallel to the leeward plane of symmetry. These results are qualitatively 
similar to Bannik and Nebbeling experiments [20], 

The solution for k=l 44 points (circumferentially) was unstable for all 
turbulence models. 

8.4 Navier-Stokes Simulation 

Computation were carried out using the NAS CRAY-2 computer, for an 
ogive-cylinder body of revolution, having a 3.5 caliber nose and a 7 caliber 
cylindrical afterbody, at a freestream Mach number of 0.2, a Reynolds 
number of 200,000 based on the body diameter, and an angle of attack of 
40°. This body geometry and test conditions correspond to those of an 
experimental study carried out by Lamont [6,7] in the NASA Ames 1 2-foot 
Pressure Wind Tunnel, in which surface pressure measurements were 
obtained. Recent analysis of the data, carried out by Hall [32], showed 
that for this laminar flow condition a wide variety of side force values 
could be obtained, depending on the orientation which the (nominally 
axisymmetric) body had in the wind tunnel. 

The computational grid (which is shown in Figure 20), consisted of 59 
points in the axial direction, 120 points circumscribing the body in the 
circumferential direction, and 50 radial points. In the computations the 
flow was assumed to be laminar; that is, no turbulence model was 
employed, and the viscosity coefficients were obtained using Sutherland’s 
viscosity law and the Stokes hypothesis. 

Three cases have been computed, all ata= 40°. In the first, the flow was 
set to free-stream conditions throughout the computational mesh, and the 
solution advanced in time. After an initial transient, the computed flow 
was found to evolve to a time-periodic, symmetric state. The second case 
was identical to the first.execpt that a very small jet was intruduced into 
the computation on one side of the body to break the symmetry of the 
solution. The jet was located axially about half length of the ogive nose, 
and was located circumferentially 90° from the angle of attack plane. In 
the third case a small geometrical disturbance (about 2% of the body 
diameter) was added to the body at the same location as the jet. In the 
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second case the solution again evolved to a time periodic state, but the 
solution was asymmetric. 

Figs. 21-23 are sample results for the symmetric computation of the 
second case. Figure 21 presents density contours in crossflow planes 
normal to axis of the body. The position of the vortex cores are indicated 
by the regions of low density, shown by the closed contours. Although the 
computational region extends circumferentially completely around the 
body, in this case the flow was found to be symmetric about the angle of 
attack plane, and thus contours for only one half of the flowfield are 
shown. One main vortex (pair) is visible about the leeward side. In 
addition, a smaller vortices are seen in the figure, originating near the 
flanks of the body on the leeward side. These small vortices are seen to 
move upwards from the body surface with increasing time, and to merge 
with the main vortex pair. As this occurs. Additional small vortices are 
observed to form at the flanks (below the upward moving vortices). This 
cycle continues, and has a definite periodic frequency. 

Time histories of normal force and pitching-moment coefficients 
(moments taken about the nose) are shown in Figs. 22 and 23, respectively. 
These show part of the initial transient, and the eventual evolution of the 
solution to a periodic stat. The period of the force and moment history 
corresponds to the time between the generation of successive small 
vortices at the flanks of the body. The occurence of the time periodic 
solution is extrimely intriguing. At present we believe that this 
phenomenon is akin to the shear-layer instability observed by Payne et al. 
[33] and Blackwelder and Gad-el-Hak [34] in experimental investigations 
of the flow about delta wings at law Reynolds numbers. As such, it may 
indicate that Navier-Stokes computations, using reasonable grid 
resolution, may permit direct computations of flow instabilities and 
initial transition to turbulence for flow about complete configurations. 
Additional numerical investigation of this phenomena is underway. 


The corresponding results for the asymmetric cases are shown in Figs. 
24-30. Figure 24 shows color contours of density in crossflow planes 
normal to the axis of the body. The position of the vortex cores are 
indicated by the regions of low density (denoted by reds and yellows). 
Dark blue contours denote densities approaching free-stream density. At 
least four major vortices are visible, which shed from the body and extend 
into the leeward-side flow. The first leaves the near side of the body at 
the ogive -cylinder junction, as indicated by the change in the vortex core 
density from red to green to blue in this region. A second vortex leaves 
the far side of the body mid-way down the cylinder, while a third leaves 
the near side of the body at the rear of the cylinder. A fourth major vortex 
exists on the far side of the cylinder, but it cannot be seen in this figure, 
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since it lies behind the near-side vortex. In addition to the four main 
vortices, small-scale vortices, similar to those observed in the symmetric 
case, are visible. 

Figure 25 shows particle traces for the same flowfield, and confirms the 
presence of multiple vortices. The particles emanating from the near side 
of the body are colored blue, while those from the far side are red. The 
vortices are seen to extend downstream over the entire body from their 
point of origin. 

Time histories of the normal force, pitching moment, and yawing moment 
coefficients (moments taken about the nose) are shown in Figs. 26-28 
respectively. As was presented for the symmetric case, these figures 
show part of the initial transient, and the evolution to a periodic state. 
The mean values of the normal force and pitching moment coefficients are 
close to those obtained for the symmetric case. However, yawing moment 
coefficient is seen to undergo a periodic variation about a non zero value, 
at the same frequency of the variation in the pitching moment coefficient. 
In contrast, in the symmetric case the yawing moment coefficient was 
zero at all times, since the solution, although unsteady, was always 
symmetric. Comparison of side forces acting on the body for the 
asymmetric case, with the net side forces caused by the jet alone, showed 
that the latter are only 6 % of total side forces. This indicates that the jet 
merely breaks the symmetry of the solution. 

Figs. 29-30 show density contours in crossflow planes normal to the axis 
of the body for the third case, where a small geometrical disturbance had 
added to the body. Although the solution has not been converged yet to a 
periodic state as in the above case, and these results are part of the 
initial transient, it is clear that the solution has been becoming 
asymmetric. We believe that this asymmetric flowfield will be converged 
to a non-symmetric periodic state similar to the one obtained when the 
symmetry of the solution was broken by a jet. 
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Fig. I - Contours of constant a, in the plane of 11 (the mass flow rate) 
vs. (dp/dr) w (the pressure gradient on the walls), where the 

pressure decreases, in the downstream direction. M c Curve - 
stability line, and zones of stability and instability. 
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Fig. 2 - Contours of constant Re, in the plane of (dp/dr) w (the pressure 

gradient on the walls) vs. <x where the pressure decreases, in the 
downstream direction. 
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Fig. 4- Zones of existence of solutions and their stability, where the 
pressure decreases, in the downstream direction and the curve 
of M c . 
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Contours of constant <x, in the plane of M (the mass flow rate) 
vs. k, where the pressure decreases, in the downstream 
direction. 
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Fig. 7 - Contours of constant <x, in the plane of M (mass flow rate) vs. 

(dp/dr) w (the pressure gradient on the walls), for the case 

where only one range of positive velocity exists, and the 
pressure increases, in the downstream direction. 



Fig. 8 - Curves of ct = 20° in the plane of M (the mass flow rate) (vs. 

(dp/dr) w (the pressure radient on the walls) for single double 

and triple positive direction velocity zones, where the pressure 
increases, in the downstream direction. 
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Fig. 9- The zones where on solution and multiple solutions exist, in 
the plane of (dp/dr) w (the pressure gradient on the walls) vs. a, 

where the pressure increases, in the downstream direction. 
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Fig. 1 2 - Conical solution. Velocity vectors, <x « 22.6° 
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Fig. 19.- Marching solution. Velocity vectors at different stations. 
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Fig. 20 - Computational grid and coordinates. 





Fig. 22 - Normal force coefficient history, symmetric flow. 







Fig. 23 - Pitching moment coefficient history, symmetric flow. 




Density contours showing vortex asymmetry. 




Part ical traces above leeward side, showing vortex asymmetry 
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Fig. 29 - Density contours. Left side of the body. Asymmetric case with 
geometrical disturbance. 



geometrical disturbance. 


